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ABSTRACT 

The twisted magnetospheres of magnetars must sustain a persistent flow of 
electron-positron plasma. The flow dynamics is controlled by the radiation field 
around the hot neutron star. The problem of plasma motion in the self-consistent 
radiation field is solved using the method of virtual beams. The plasma and radi- 
ation exchange momentum via resonant scattering and self-organize into the "ra- 
diatively locked" outflow with a well-defined, decreasing Lorentz factor. There is 
an extended zone around the magnetar where the plasma flow is ultra-relativistic; 
its Lorentz factor is self-regulated so that it can marginally scatter thermal pho- 
tons. The flow becomes slow and opaque in an outer equatorial zone, where the 
decelerated plasma accumulates and annihilates; this region serves as a reflector 
for the thermal photons emitted by the neutron star. The e ± flow carries electric 
current, which is sustained by a moderate induced electric field. The electric 
field maintains a separation between the electron and positron velocities, against 
the will of the radiation field. The two-stream instability is then inevitable, and 
the induced turbulence can generate low-frequency emission. In particular, ra- 
dio emission may escape around the magnetic dipole axis of the star. Most of 
the flow energy is converted to hard X-ray emission, which is examined in the 
accompanying paper. 

Subject headings: plasmas — stars: magnetic fields, neutron — X-rays 

1. Introduction 

The observed activity of magnetars is believed to be caused by their surface motions, 
which are driven by strong internal stresses. The magnetosphere is anchored in the neutron 
star and twisted by the surface motions, resembling the behavior of the solar corona. As a 
result it becomes non-potential, V x B / 0, and threaded by electric currents (Thompson 
et al. 2002; Beloborodov & Thompson 2007). The currents flow along B, i.e., the twisted 
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magnetosphere remains nearly force-free, j x B = 0. First quantitative models of dynamic 
twisted magnetospheres of neutron stars have been calculated recently by Parfrey et al. 
(2012a,b). These calculations demonstrate how the twist creates spindown anomalies and 
initiates giant flares when the magnetosphere is "overtwisted" and loses equilibrium. 

The free energy stored in a magnetospheric twist is gradually dissipated through the 
continual electron-positron discharge that sustains the electric current. As a result, the 
magnetosphere tends to "untwist" on the characteristic ohmic timescale of a few years. 
Electrodynamics of untwisting is quite peculiar: whenever the crustal motions stop or slow 
down, the electric currents tend to be quickly removed from the magnetospheric field lines 
with small apex radii -R max (Beloborodov 2009). Currents have longest lifetimes on field 
lines with -R max ^> R (where R ~ 10 km is the star radius) and form an extended "j- 
bundle" (Figure 1). The j-bundle has a sharp boundary, which gradually shrinks toward 
the magnetic dipole axis. In particular, the footprint of the j-bundle on the neutron star 
surface, which may be observed as a hot spot, shrinks with time. Shrinking hot spots were 
indeed reported in "transient magnetars" whose magnetospheres were temporarily activated 
and then gradually relaxed back to the quiescent state (see Figure 5 in Beloborodov 2011a 
and refs. therein). 

The j-bundle must be filled with relativistic plasma that carries the electric current 
j = (47r/c)V x B. The plasma is continually created by discharge near the star and must 
expand along the extended field lines. The plasma emits persistent nonthermal emission, 
converting the dissipated twist energy to radiation. In the accompanying paper (Beloborodov 
2012, hereafter B12), we calculate the spectrum of produced radiation and show that it forms 
the observed hard X-ray component with a peak around 1 MeV. The present paper studies 
in more detail the dynamics of the e ± plasma. 

Previously, semi-transparent plasma flows around magnetars were invoked to explain 
the deviation of the observed 1-10 keV emission from a thermal spectrum (Thompson et 
al. 2002). It was usually assumed that the magnetar corona is filled with positive and 
negative charges that are counter-streaming with mildly relativistic speeds. The counter- 
streaming picture was motivated by the fact that electric current j must flow along the 
twisted magnetic field lines. The coronal plasma must be nearly neutral; it can carry the 
required current if the opposite charges with densities n + = n_ flow in the opposite directions, 
so that j = e(v + n + — where v + V- < 0. The velocities v± are free parameters in this 

phenomenological model, which may be adjusted so that resonant scattering of thermal X- 
rays in the corona reproduces the 1-10 keV part of the magnetar spectrum (e.g. Fernandez 
& Thompson 2007; Nobili et al. 2008; Rea et al. 2008). 

The counter-streaming model has, however, a problem. Note that the thermal radiation 
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Fig. 1. — Snapshot of a slowly untwisting magnetosphere. In this example, a global twist 
with a uniform amplitude ip = 0.2 was implanted into the dipole magnetosphere at t — 0, 
and the snapshot shows the magnetosphere at t ~ 1 yr. Details of the calculations are 
described in Beloborodov (2009). The plane of the figure is the poloidal cross section of the 
magnetosphere. The black curves are the poloidal magnetic field lines. The magnetosphere is 
symmetric about the vertical axis and the equatorial plane; therefore, the figure only shows 
one quarter of the poloidal cross section. The neutron star is shown by the black circle 
(radius R ~ 10 km). Left panel: Current density j normalized to BR/c. The region from 
which currents have been pulled into the star (the potential "cavity" with j = 0) is shown 
in white. The boundary between the cavity and the j-bundle (magenta curve) expands with 
time, i.e. the j-bundle shrinks toward the vertical axis. Right panel: Twist amplitude ip at 
the same time. The twist amplitude is defined for each closed field line as the azimuthal 
displacement of its footpoint in the southern hemisphere relative to the footpoint in the 
northern hemisphere. 
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of the magnetar [hv ~ 1 keV) is resonantly scattered at large radii r ~ 10R fa 100 km where 
B ~ 10 11 — 10 12 G. 1 In this region, the plasma is strongly pushed by radiation away from the 
star and the counter-streaming model needs an electric field E» that forces charges of the 
right sign to move toward the star against the radiative drag. At the same time, the electric 
field acts on the opposite charges and accelerates them away from the star, cooperating 
with the radiative push. In the presence of plasma (which is inevitably created near 
magnetars), the outward acceleration generates relativistic particles, and no self-consistent 
solution exists for the mildly-relativistic counter-streaming model. 

In this paper (and the accompanying paper B12), we develop a different picture of 
plasma circulation in the magnetar corona. It is schematically shown in Figure 2. The outer 
corona is inevitably filled with e ± pair plasma of a high density n, which is larger than j/ec 
by the "multiplicity factor" M 3> 1; in this respect it resembles the flow along the open 
field lines of a rapidly rotating, strongly magnetized neutron star (e.g. Hibschman & Arons 
2001; Thompson 2008a; Medin & Lai 2010). Both electrons and positrons outflow from 
the magnetar, and radiation pressure forces the particles to accumulate in the equatorial 
plane of the magnetic dipole, where they annihilate. The required current j = (c/47r)V x B 
is sustained in the outflow by a moderate electric field En. This field is self-consistently 
generated to maintain a small difference between the velocities of the ± charges, (v+ — 
v-)/v± ~ M.^ 1 <c 1, so that the condition e{n + v + — ri-vJ) = j is satisfied with n + ~ n_ 
and v + V- > 0. In the simplest, two-fluid model (Section 3) the velocities v± tend to be 
"locked" by the balance of two forces, electric and radiative. 

The coronal outflow significantly changes the radiation it interacts with via scattering. 
The problem of outflow dynamics can be formulated as a problem of self-consistent radiative 
transfer where particles and photons exchange energy and momentum as they flow away from 
the neutron star. This problem is solved in this paper using a specially designed numerical 
method. 

The paper is organized as follows. Section 2 discusses the creation of e ± pairs and 
their circulation in the inner and outer magnetosphere. Section 3 presents the model of a 
radiatively locked outflow in its simplest version using a two-fluid description and assuming 
an optically thin magnetosphere. Section 4 discusses the two-stream instability in the 
flow and the origin of low-frequency emission from magnetars. Then, in Sections 5 and 6 
we formulate and solve the full problem where an outflow with a broad momentum distri- 
bution function and significant optical depth interacts with the neutron-star radiation. The 
numerical method and results are described in Section 6. Our conclusions are summarized 



This assumes that photons are scattered by electrons, not ions. 
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in Section 7. 



2. Creation and circulation of e ± 



pairs 



The creation of e ± pairs by an accelerated particle is a two-step process: the particle 
generates a high-energy photon (via resonance scattering) and then the photon converts to 
an e ± pair in the strong magnetic field. An accelerated electron (or positron) can resonantly 
scatter photons of energy fkj once it reaches the Lorentz factor required by the resonance 
condition 7(1 — (3cos$)oj = ug, where 1? is the angle between the photon and the electron 
velocity and ojb = eB/m e c. When 1? is not small, the resonance condition gives 



The scattered photons are boosted in energy by the factor of ~ 7g C . Such high-energy 
photons quickly convert to e ± pairs in the strong magnetic field, creating more particles near 
the star. A similar process of e ± creation operates in the polar-cap discharge of ordinary 
pulsars, but in a different mode. In ordinary pulsars, the high-energy photons convert to 
e ± with a significant delay. The scattered photon initially moves nearly parallel to B and 
converts to e ± only when it propagates a sufficient distance where its angle # 7 with respect 
to B increases so that the threshold condition for conversion is satisfied. This delay leads 
to the large unscreened voltage in pulsar models. 

In contrast, the magnetic field of magnetars is so strong that pair creation can occur 
immediately following resonant scattering (Beloborodov & Thompson 2007). The energy of 
a resonantly scattered photon is related to its emission angle # 7 by 



where E B = (2B/B Q + l) 1 / 2 m e c 2 is the energy of the first excited Landau level and Bq = 
mlc 3 /he ~ 4.4 x 10 13 G. The scattered photon may immediately be above the threshold for 
conversion, E > E t ^ r = 2m e c 2 / sin# 7 , if B > 4Bq. Therefore, discharge in magnetars can 
screen E\\ more efficiently than in ordinary pulsars and buffer the voltage growth once the 
Lorentz factors of accelerated particles reach 7 SC given in Equation (1). 
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2.1. Pair creation on field lines with apexes i? max ^ 2i? 

The discharge on twisted closed field lines can be explored using a direct numerical 
experiment where plasma is represented by a large number of individual particles in the 
self-consistent electric field. The existing numerical simulations (Beloborodov & Thompson 
2007) describe the discharge on field lines that extend to a moderate radius R max ^ 2R, 
where R is the radius of the neutron star. The magnetic field is ultrastrong everywhere 
along such field lines, B 3> Bq, and resonant scattering events may effectively be treated as 
events of pair creation — a significant fraction of scattered photons immediately convert to 
e ± . 

The simulations demonstrate that voltage and pair creation self-organize in the twisted 
magnetosphere so that a particle on average scatters ~ 1 photon as it travels through the 
electric circuit, maintaining the near-critical multiplicity of pair creation A4 ~ 1. This 
criticality condition regulates the induced voltage to $ e ~ 10 9 V, which accelerates 
particles to Lorentz factors 7 ~ 10 3 . The electric circuit operates as a global discharge, in 
the sense that the accelerating voltage is distributed along the entire field line between its 
footpoints on the star. It is quite different from the localized "gap" that is usually considered 
above polar caps in pulsars. 

The discharge fluctuates on the light-crossing timescale ~ Rj c and persists in the state of 
self-organized criticality. The behavior of the global circuit resembles a continually repeating 
lightning: voltage between the footpoints of the field line quasi-periodically builds up and 
discharges through the enhanced production of charges. The average plasma density in the 
circuit n is close to the minimum density n min = j/ec, as required by the criticality condition 
M = n/n min ~ 1. 

The global-discharge picture applies only to field lines with -R max ~ 21? and becomes 
irrelevant when the currents are erased in the inner magnetosphere as shown in Figure 1. 
Then the observed activity must be associated with currents on field lines with large -R max , 
i.e. extending far from the star. 

2.2. Pair creation on field lines with apexes R mSLX 3> R 

The discharge on extended field lines may be expected to have a similar threshold 
voltage $ e ~ 10 9 V, because the conversion of upscattered photons to is efficient near 
the field-line footpoints where B 3> Bq. In this zone, particles are able to resonantly scatter 
soft X-rays once they are accelerated to 7 ~ 10 3 (Equation 1), which requires $ e ~ 10 9 V. 
Further growth of voltage would cause excessive creation of e ± moving in both directions, 
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toward and away from the star, leading to efficient screening of E\\. 

A large fraction of the created particles must outflow tor > along the extended field 
lines. The rate of resonant scattering by a relativistic particle increases as it moves from 
B 3> Bq to B ^ Bq. The particle scatters many more photons, because the resonance 
condition shifts toward photons of lower energy hcu ies oc B whose number density is larger. 
Note also that the effective cross section for resonant scattering a rcs = 2ir 2 r e c/uj res increases 
as B~ x (here r e = e 2 /m e c 2 is the classical electron radius). Practically all photons scattered 
by the outflowing particles in the region B ^ 10 13 G convert to e ± ; detailed Monte-Carlo 
simulations of this process are presented in the accompanying paper (B12). In essence, 
the particles outflowing from the discharge zone lose energy to photon scattering, and this 
energy is transformed to new generations of e ± . As a result, the e ± multiplicity of the outflow 
increases from M ~ 1 to M ~ 100. This implies that there is no charge starvation in the 
outer corona — there are plenty of charges to conduct the current demanded by the twisted 
magnetic field. 

Pair creation sharply ends near the surface of B ps 10 13 G; outside this surface the 
resonantly scattered photons are not absorbed. The steady relativistic outflow without pair 
creation maintains M. = const along the magnetic field lines. This follows from conservation 
of magnetic flux, charge, and particle number, which give n/B = const, j/B = const, and 
n/j — const along the field line. 

2.3. Global circulation of pair plasma 

The picture of plasma circulation in the magnetar magnetosphere is summarized in 
Figure 2. The inner corona (field lines with i? max of a few magnetar radii) is filled with 
the ultra- relativistic counter-streaming e~ and e + . In this region, resonant scattering is 
marginally efficient and a global discharge operates as described in Section 2.1, with multi- 
plicity A4 ~ 1. Outside this region, pair multiplicity is much higher and the electric current is 
organized with both e~ and e + outflowing from the star. The exact location of the boundary 
between the two regions depends on the strength of the magnetic field of the star. 

In the outer corona, the opposite flows in the nothern and southern hemispheres meet in 
the equatorial plane of the magnetic dipole and stop there. Two effects prevent their inter- 
penetration. (1) Radiative drag is strong in the outer corona and pushes both nothern and 
southern flows toward the equatorial plane (see Section 3 below). (2) When the two opposite 
flows try to penetrate each other, a two-stream instability develops. As a result, a strong 
Langmuir turbulence is generated, which inhibits the penetration. This effect is particularly 
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Fig. 2. — Schematic picture of plasma circulation in the magnetosphere with surface B ~ 
10 15 G. Two regions are indicated. (1) "Inner corona." Here e have a moderate multiplicity 
Ai ~ 1. The particles do not stop in the equatorial plane. The electric field E» ensures that 
electrons and positrons circulate in the opposite directions along the magnetic field lines, 
maintaining the electric current demanded by V x B. The particles are lost as they reach 
the footpoints of the field line and continually replenished by pair creation. (2) "Outer 
corona" — extended field lines with R max 3> R. Electrons and positrons are created by the 
discharge near the star and some of them flow outward to the region of weaker B. Here 
resonant scattering enhances the pair multiplicity, M. ^> 1, and decelerates the outflow. The 
e ± particles stop at the apexes of magnetic field lines (blue region in the equatorial plane), 
accumulate, and annihilate there. The number fluxes of electrons and positrons toward the 
annihilation region differ by a small fraction ~ .M -1 , so that the outflow carries the required 
electric current j = (c/Att)V x B. Electrodynamics of the twist dissipation implies that the 
inner corona is less likely to be active, as the electric currents are erased by the expanding 
cavity (Figure 1); the observed activity tends to concentrate on extended field lines that 
form the outer corona. 
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important in the transition region between the outer and inner corona. For these field lines, 
resonant scattering is efficient enough to generate M. > 1 but not strong enough to stop 
the pair plasma in the equatorial plane. Then the colliding northern and southern flows are 
stopped by the two-stream instability. This behavior contrasts with the inner corona where 
the induced electric field enforces the counter-streaming of e~ and e + , i.e. the opposite flows 
with Ai ~ 1 are forced to penetrate each other despite the two-stream instability. 

The density of e 1 * 1 pairs accumulated in the outer equatorial region (shown in blue in 
Fig. 2) is regulated by the annihilation balance. In a steady state, the annihilation rate is 
given by iV arm 2Ai(I/e), where / is the electric current through the annihilation region. 
The corresponding annihilation luminosity is L ann = 2m e c 2 A4 I/e. 



3. Radiatively locked coronal flow 

Dynamics of the e ± flow in the outer corona (r ^> R) is influenced by resonant scattering, 
which exerts a strong force J 7 on the particles along the magnetic field lines. J 7 vanishes only 
if the particle has the "saturation momentum" p± such that the radiation flux measured in 
the rest frame of the particle is perpendicular to B. In the simplest case of a weakly twisted 
dipole magnetosphere exposed to central radiation p± is given by (Appendix B), 2 

. „. 2 cos 6* 

p.M) = — T , (3) 

where momentum is in units of m e c. The radiative force always pushes the particle toward 
p = p±. The strength of this effect may be measured by the dimensionless "drag coefficient," 

rT 

V = -— r (4) 

pm e c z 

Momentum p+ is a strong attractor in the sense that deviations p — p+ generate V ^> 1 in 
the outer corona (see Appendix A). 

Even an extremely strong radiative drag does not imply that e + and e~ acquire exactly 
equal velocities (5 + = (5- = f3+ = p^l+pD^ 1 / 2 . Such a "single-fluid" flow would be unable to 
carry the required electric current j, and E\\ must be induced to ensure a sufficient velocity 
separation f3 + — /?_. Below we describe the two-fluid model with /3 + 7^ /3_. 



2 This expression is valid in the region where 1 — B r /B > (R/r) 2 . In this region, stellar radiation may be 
approximated as a central flow of photons, neglecting the angular size of the star ~ R/r. 
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3.1. Two-fluid model: basic equations 

Consider an e^ 1 flow in the region where no new pairs are produced. In a steady state, 
the fluxes of electrons and positrons are conserved V • (n± v±) = 0. This implies V • j ± = 
where j ± = rken±v± are the contributions of electrons and positrons to the net electric 
current j = j + + j . Since e ± move along the magnetic field lines, j ± = a± B where a± are 
scalar functions. From V • j ± = and V • B = one gets B • Va± = 0. Therefore, 

^- = const, = const, (5) 

are constant along the field line. The net current j = j + + j_ is fixed by the condition 
j = (c/47r)|V x B|. The multiplicity of pairs is defined by 

M = j -±±\t\. (6) 
J 

Here "+" corresponds to positrons, which carry current j + > and "— " corresponds to 
electrons, which carry j_ < 0; we assume a positive net current j = j + + j_ for definiteness. 
A counter-streaming model (v + v_ < 0) would have At — 1. A charge-separated outflow 
(n_ = 0) would also have M. — 1. We consider here pair- loaded outflows with M. > 1. 

The outflowing e ± plasma must be nearly neutral, 3 

n + ?s n_, (7) 

otherwise a huge electric field would be generated that would restore neutrality. Using this 
condition, one finds that en + v + — en_v_ = j is satisfied if 

A deviation from this condition implies a mismatch between the conduction current j + + 
j_ and (c/47r)V x B, which would induce a growing electric field according to Maxwell 
equation d~E/dt = cV x B — An j. An electric field E\\ must be established in the outflow to 



3 Here we neglect rotation of the neutron star and its magnctosphere, which is a good approximation 
everywhere except the open field-line bundle that connects the star to the light cylinder. In a more exact 
model, the Gauss law in the co- rotating frame V • E = An(p + p v ) includes the effective vacuum charge 
density p v = £1 ■ H/2ttc where f2 is the angular velocity of the star (Goldreich & Julian 1969). Then the 
neutrality condition becomes e(n + — n_) + p v =0. Magnetars rotate slowly (typical O <~ 1 rad/s) and hence 
a small fraction of their magnetic flux is open, typically ^ 1CT 4 . In the main, closed magnetosphere, the 
condition \p v /e\ -C \ j/ec\ < n± is satisfied, and neutrality requires n + w n_ with a high accuracy. 
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sustain the condition (8) against the radiative drag that tends to equalize (3- and (3 + at f3+. 
This moderate electric field must be self-consistently generated by a small deviation from 
neutrality, 5n = n + — n_ <C n± . 

The two-fluid dynamics of the outflow is governed by equations, 

where / is length measured along the magnetic field line. In the region of strong drag, 
\T>\ 3> 1, the left-hand side is small compared with J 7 ; then the radiative and electric forces 
on the right-hand side nearly balance each other, -F(7+) ~ — eE\\ and ^(j-) ~ eE\\. This 
implies, 

^( 7+ ) ~ ~Hl-) {\V\ » I)- (10) 
Equations (8) and (10) describe the "radiatively locked" two-fluid current with pair multi- 
plicity M.. The solution 7 ± to these two equations exists if M. > 1. Note that j3 + 
in the radiatively-locked state. 

Let us now relax the assumption \T>\ ^> 1. Then the inertial term m e c 2 d^±/ dl must 
be retained in Equations (9), i.e. we now deal with differential equations for 7 ± . Since 7_ 
and 7+ are not independent — they are related by condition (8) — it is sufficient to solve 
one dynamic equation, e.g. for 7+ (and use the dynamic equation for 7_ to exclude £j|). 
Straightforward algebra gives, 

2 di+ ^(7+) + F(n-) (M - 1 \ 2 (l- 



meC dl ~ l + rf 7 -/rf7+ ' d 1+ ~\M + l) ■ 



3.2. Sample numerical model 

Suppose that e ± plasma is injected near the star with a given multiplicity, e.g. M. = 50, 
and a given high Lorentz factor, e.g. 7+ = 100. The corresponding 7_ is determined by 
Equation (8). Suppose that the plasma is illuminated by the blackbody radiation of the 
star of temperature kT = 0.5 keV and neglect radiation from the magnetosphere itself. This 
approximation is valid only for optically thin magnetospheres, which are considered here for 
simplicity. A more detailed model will be developed in Sections 5 and 6, which will take into 
account radiation scattered in the magnetosphere; then radiation field is not central, and we 
will have to solve radiative transfer to determine the flow momentum. 

An explicit expression for ^(7) exerted by the central radiation is given in Appendix B 
(Equation B6). The steady-state solution for 7± in the outer corona can be found by inte- 
grating Equation (11) along the magnetic field lines. The result is shown in Figure 3. The 
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Fig. 3. — Lorentz factors 7+ (left panel) and 7_ (right panel) in the two-fluid model of the 
flow. In this example, the electric current is carried by the e outflow of a fixed multiplicity 
M. = 50. The plasma is injected at radius r = 2R and outflows along the magnetic field 
lines (white curves). The flow is illuminated by the star with temperature kT = 0.5 keV 
(magenta circle at the origin), and the radiation exerts the forces J r {'j±) on the positron 
(+) and electron (— ) fluids. The Lorentz factors 7+ and 7_ change as the flow enters the 
drag-dominated region \T>\ 3> 1. The region \T>\ > 3 is shown by the thick black curve and 
shadowed in black. V < for positrons (7+ > 7*) and V > for electrons (7. < 7*). 
The radiative drag stops the plasma in the equatorial plane outside ~ 8R. A nearly dipole 
magnetic field (weakly twisted) with -B po i e = 10 15 G is assumed in this example. R is the 
neutron-star radius. 
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relativistic outflow is injected near the star and initially weakly interacts with the radia- 
tion; then it enters the drag-dominated region \T>\ ^> 1. The solution is not sensitive to 
the precise radius of e ± injection as long as it is small enough, before the plasma enters the 
drag-dominated region. 

The electric field in the region of \V\ 3> 1 is given by eE\\ « and the corre- 

sponding longitudinal voltage established in the outer corona is found by integrating -F(7+) 
along the field line, e$ e — J -F(7+) dl. Its typical value for the model in Figure 3 is 
~ 10 7 V. Flows with lower Ai develop stronger electric fields, however in all cases of interest 
{M. 1) the drag-induced voltage is below 10 9 V. 

The calculations shown in Figure 3 assume that the magnetospheric plasma is every- 
where optically thin. This is not so for real magnetars. Thompson et al. (2002) showed 
that the characteristic optical depth r of a strongly twisted magnetosphere (twist amplitude 
tp ~ 1) is comparable to unity. When the large pair multiplicity M. is taken into account, 
the estimate changes to 

r ~ —Z » 1. (12) 

P± 

This estimate describes the optical depth seen by photons that can be resonantly scattered 
by the flow, i.e. the resonance condition 7(1 — (3cosi9)u) = ub is saisfied somewhere along 
the photon trajectory. The large r implies the presence of scattered radiation in the mag- 
netosphere, which is quasi-isotropic rather than central. This increases the drag exerted on 
the outflow and reduces p*. One may also expect a self-shielding effect: the drag force T 
experienced by an electron (or positron) is reduced by the factor of r _1 . The problem of 
self-consistent outflow dynamics will be solved in Sections 5 and 6. The main feature seen 
in Figure 3 will persist in the full self-consistent solution: the outflow is strongly decelerated 
(and drag-dominated) in the equatorial region at r ~ (10 — 20) i?. 



4. Two-stream instability, anomalous resistivity, and radio emission 

4.1. Two-stream instability 

The two-fluid flow with v + > t>_ is prone to the two-stream instability (e.g. Krall & 
Trivelpiece 1973). The growth rate of the instability is obtained from the dispersion relation 
for Langmuir modes with frequency u and wavevector k, which can be derived by considering 
perturbations 5 r y± of the two-fluid system and using continuity, Euler, and Poisson equations, 

U) 2 , U) 2 

1 -0, (13) 



7^(w — kv + ) 2 — kvJ)'' 
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where u± = 47rn±e 2 /m e . The plasma is nearly neutral, n + = and hence u + = oj-. 
Equation (13) has a solution uj(k) with a positive imaginary part that describes an unstable 
mode. The solution simplifies in the following two limits: 

(1) 7+ — 7_ <C 7±, which is valid when M 3> 7± (see Equation 8). Then the flow is 
convenient to view in its center-of-momentum frame that moves with 7 « (7+ + 7-)/2. In 
this frame, the two fluids with densities n± = 7 _1 n± move in the opposite directions with 
non-relativistic velocities v± = ±v, so the dispersion relation (13) simplifies. It gives the 
most unstable mode k = (y/3/2) oo±/v with a growth rate Y = u±/2. 

(2) 7+/7_ 1. Then the contribution of positrons to the dispersion relation is small 
compared to that of the electrons. The growth rate of the instability is given by T 
7_ 7+ w p (e.g. Lyubarsky & Petrova 2000). 

This estimate gives the characteristic length-scale of the instability c/Y that is much 
shorter than the electron free path to resonant scattering, A sc . Hence the radiation drag and 
the induced electric field eE\\ = ±J r (p±) are unable to lock the positive and negative charges 
at the momenta p + and p_ calculated in the two-fluid model. The instability will generate 
plasma oscillations that should broaden the momentum distribution so that particles fill the 
region p_ < p < p + . 

The generated plasma oscillations may be expected to introduce an anomalous resisti- 
tivity. The fluctuating En in the oscillations creates a stochastic force that tends to reduce 
the free-path of a charged particle. A simplest estimate suggests that this effect could be 
very strong. Suppose a substantial fraction of the energy density of the flow / ym e c 2 n is 
given to plasma oscillations. Then the characteristic electric field is En ~ (87T7m e c 2 n) 1 / 2 ; it 
is much stronger than En in the radiatively locked two-fluid model, however it is irregular 
and quickly changes sign. Suppose that the stochastic electric force exerted on the particle 
randomly changes sign on a timescale At ~ ^p 1 , where u p is the plasma frequency. The 
stochastic En gives the momentum kicks AP ~ eE\\u~ l and causes diffusion of particles in 
the momentum space with the diffusion coefficient 

n (AP)2 (eE|l)2 (u) 

Dp ~ - a* — ^r- (14) 

Diffusion in momentum space p 2 (t) ~ D p t implies a small free-path of the particle, A ~ 
(7/3) 2 mgC 3 cu p /(ei?||) 2 , much smaller than the mean free path to resonant scattering. Thus, 
a large anomalous resistivity could, in principle, be possible, and then a large longitudinal 
voltage would be generated to maintain the electric current. 
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4.2. Numerical experiment 

To explore the role of the two-stream instability and anomalous resistivity, we designed 
the following numerical experiment. Keeping in mind that particles around magnetars can 
flow only along the magnetic field lines, consider the simple one-dimensional problem. Sup- 
pose an beam is continually injected at the boundary z = of the computational box 
< z < L. The rate of electron injection is smaller than the rate of positron injection, so 
that the flow carries current j > 0. Positrons are injected with fixed p + and electrons with 
fixed p_ < p + ; we chose p + = 7 and p_ = 0.5. The simulation keeps track of the flow in the 
computational box of length L ~ 150c/u p . The escape boundary condition is implemented 
at z — L. The flow is simulated as a large collection of individual particles (N ~ 10 6 ) that 
move in their collective electric field (see Beloborodov & Thompson [2007] for details of the 
numerical method). The simulation was run for a time of t ~ lOOL/c, long enough to see 
the quasi-steady behavior. 

The results of the simulation are as follows. As expected, strong plasma oscillations 
develop and persist in the flow, as the instability is continually fed by the injection of the 
two streams at z — 0. The fluctuating electric field reaches very high amplitudes that could 
stop a particle with Lorentz factor 7+ on a short scale, much shorter than L. A snapshot of 
En (z) is shown in Figure 4. The integral of the electric field over z determines the voltage $ e 
between the two boundaries of the computational box. The measured <3> e in the simulation 
fluctuates in time (Figure 5), and we calculated its value averaged over lOOL/c. This value 
turns out to be very small, e$ e ~ 0.l7 + m e c 2 . Thus, the measured anomalous resistivity is 
small, in sharp contrast with the simplest estimate (14) that would predict A <C L and hence 
e$ e ^> 7 + m e c 2 . 

The failure of the estimate (14) is related to the assumption that the stochastic electric 
force applied to the particle is random, uncorrelated on timescales longer than oo^ 1 . The 
numerical simulation indicates that this assumption is incorrect. Apparently, a complicated 
time- dependent pattern is organized in the phase space, which allows the charges to find 
small-resistance paths through the waves of E\\ and conduct the current at a low net voltage 
and a low dissipation rate. 

One limitation of the numerical experiment should be noted: the one-dimensional com- 
putational box allows no dependence of En on the transverse coordinates x and y, which 
excludes the coupling of Langmuir waves to the transverse electro-magnetic modes. A two- 
dimensional (or a complete three-dimensional) simulation will be needed to explore the role 
of the transverse modes. We anticipate that a low anomalous resistivity will be found in 
the full simulations. A high resistivity would imply a quick dissipation of magnetospheric 
currents, which would produce a high luminosity and quickly erase the magnetic twist. This 
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Fig. 4. — Snapshot of electric field in the simulation box. 
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Fig. 5. — Voltage across the simulation box as a function of time. 
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is not supported by observations of magnetars. The observed luminosities and evolution 
timescales are consistent with the model neglecting the anomalous resistivity, where voltage 
is controlled by the threshold of the discharge, e$ e ~ 10 9 V. 

The numerical simulation shows that the two-stream instability significantly changes 
the momentum distributions of electrons and positrons from the injected delta-functions 
5(p — p ) and S(p — p + ). The e ± distributions are broadened so that they fill the region 
between the injection momenta p + and p_ (Figure 6). 

4.3. Low-frequency emission 

An important implication of the two-stream instability is the excitation of a strong 
plasma turbulence that can generate coherent low-frequency radiation. The two-stream 
instability is often considered in pulsar models as a mechanism feeding radio emission from 
the open field- line bundle (e.g. Sturrock 1971; Cheng & Ruderman 1977). A related model 
invoking radiative drag was considered by Lyubarsky & Petrova (2000). 4 In the flows 
around magnetars, the two stream instability is naturally driven by the strong electric current 
in the system that tends to lock itself in the two-fluid configuration as described in Section 3. 
The instability is continually pumped by the radiative drag in the dense radiation field of 
the magnetar. The generated low-frequency radiation has the best chance to escape near 
the magnetic axis, where the plasma density is lowest and its Lorentz factor is highest. Note 
that the two-stream instability operates on the closed (twisted) field lines, which allows the 
low-frequency emission to be bright, even though the voltage $ e ~ 10 9 V is small by the 
ordinary-pulsar standards. 

Radio pulsations have been detected and studied in detail in two magnetars XTE 1810- 
197 and IE 1547.0-5408 (Camilo et al. 2006; 2007). The estimated radio luminosity L r ~ 
10 30 erg s _1 requires a sufficiently high ohmic dissipation rate, IQ e > L r , which cannot be 
generated in the open field-line bundle unless $ e exceeds L r //Lc ~ 10 n £ r ,3o V. Here /lc is the 
electric current circulating in the open bundle, Jlc ~ c /V-^lc> wnere -^lc — cP/2-k ~ 10 10 cm 
is the light cylinder for a magnetar rotating with the period P of a few seconds. As we argued 
above, voltage is likely near the threshold of e ± discharge, which gives $ e <C 10 11 V. Then 



4 They suggested that a broad momentum distribution of relativistic particles in the open field-line bundle 
can evolve into a two-hump distribution as a result of resonant scattering losses, as the lower energy particles 
lose their energy faster than the more energetic ones. In contrast, the two-fluid flow described in Section 3 is 
shaped by the induced E\\ so that it sustains the electric current j; without E\\ radiative drag would equalize 
the velocities of all particles at v*. 
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the electric current associated with observed radio emission, I ^> Jlc, must flow in the 
closed (twisted) magnetosphere and give a bright and relatively broad radio beam. This is 
consistent with the unusually broad radio pulses of magnetars, much broader than the typical 
pulse of ordinary pulsars with similar periods (Camilo et al. 2006; 2007). Note also that 
the plasma density in the twisted closed magnetosphere reaches much higher values than in 
the open field-line bundle, and the plasma frequency may approach the infrared band. This 
may explain the observed hard radio spectra. 

One could consider the possibility that the radio luminosity of magnetars is generated 
by enhanced dissipation in the open field-line bundle and its immediate vicinity. Thompson 
(2008b) suggested that the diffusion of magnetic twist ip in the closed magnetosphere initiates 
a strong Alfvenic turbulence near the light cylinder, with a high dissipation rate. This picture 
assumes that the magnetic twist tends to spread due to ohmic diffusion, as observed in normal 
laboratory plasma with a finite resistivity. However, later work (Beloborodov 2009; 2011a) 
showed that the twists in neutron-star magnetospheres evolve differently: the twist is erased 
"inside out" rather than spreads diffusively. The twist evolution dip/dt is controlled by 
voltage induced along the magnetic field lines, $ e , so that dip/dt oc d^e/dJ 7 , where J 7 is 
the magnetic flux function labeling the field lines (J 7 = on the magnetic dipole axis). An 
increased voltage near the axis would imply a large negative d^e/dJ 7 , which implies a large 
negative dip/dt, i.e. rapid untwisting. Thus, the high-$ e twist near the open field-line bundle 
cannot be sustained. The ohmic effects in the closed magnetosphere can pump the twist near 
the open bundle only if <3> e is reduced toward the magnetic dipole axis, i.e. d^e/dJ 7 > 0. 
Then the twist pumping continues until the expanding cavity j = reaches the open bundle; 
a snapshot of this evolution is shown in Figure 1. The pumping of ip leads to outbursts and 
spindown anomalies (Parfrey et al. 2012a,b). 

Our preferred scenario for the low-frequency emission from magnetars may be summa- 
rized as follows. The emission is generated by the two-stream instability on the twisted closed 
field lines with the apex radii R max such that R <C R ma , x *C -Rlc- These field lines carry 
a large electric current / ~ ipcfi/R^ n&x ^> Jlc and a modest voltage $ e ~ 10 9 V. The high 
plasma density and the broad beam of radiation expected on these field lines may explain 
the unusual radio pulsations of magnetars. 



- 21 - 



5. Dynamics of outflow with a broad momentum distribution 

5.1. Waterbag model 

The plasma instability discussed in Section 4 is generated by the gradient of the dis- 
tribution function df e /dp, and the feedback of the excited plasma waves tends to make the 
distribution flatter. The numerical simulation in Section 4 illustrates how two interpenetrat- 
ing cold fluids of e + and e~ with f e {p) = (1/2) [S(p — pJ) + 5(p — p + )} quickly evolve into a 
state with a broad and smooth f e (p). Below we design a simple modification of the two-fluid 
model that takes this effect into account. 

The simplest model has a top-hat distribution function. In plasma physics, this approx- 
imation is often called "waterbag" model. Like the two-fluid model, the outflow is described 
by two parameters p± (or /3±). However, now instead of two delta-functions we require the 
e ± distribution to be flat between p_ and p + , 

feip) = { {P+ P ~ <P<P + (15) 

J yF> \ p<p„ or p>p+ y ' 

This distribution includes both electrons and positrons. The plasma must be nearly neutral, 
n + = n_ , and the total density of particles n = n + + n_ is given by 

n — -= - = — =-, 16 

13c eP' K ' 

where N = Aij/e is the particle number flux and (3 is the average velocity, 

P= f P(p)fe(p)dp=^^-. (17) 
J p+ — P- 



Similar to the two- fluid model, p + and p_ are not independent, as the outflow must carry 
current j with a given multiplicity M.; this would be impossible if, for instance, p + — p_. 
The minimum possible width of the distribution function (i.e. minimum p + — p_) is achieved 
if all negative charges are slower than all positive discharges, as shown in Figure 7. We 
will adopt this idealized momentum distribution in our numerical simulations. This is a 
rather crude approximation to more realistic distributions of e + and e~, which overlap in 
momentum space (cf. Figure 6). 5 Nevertheless, it is useful as it allows one to explore all 
basic features of the e ± outflow using a concrete relation between p + and p_. This relation 



5 Thc minimum-width waterbag model may be particularly crude near the star where the outflow just 
begins to experience significant radiative drag. The faster positive charges will experience drag first while 




Fig. 7. — Waterbag distribution of minimum width. 
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is determined by the outflow multiplicity M.. It is easy to show that the generalization of 
the two-fluid Equation (8) to e ± flows with any distribution function is given by 

where f3 + and /3„ are the average velocities of the positive and negative charges, respectively. 
For the waterbag distribution shown in Figure 7 this condition gives 

gz _ 7(P)~7- =1 2_ fig) 

^+ 7+"7(p) -M + l' ^ 2 ' liyj 

where 7(2?) = (1 + p 2 ) 1 ^ 2 . Equation (19) determines the relation between p + and p_ (for a 
given it is shown in Figure 8. 

Taking into account the relation between p + and p_, the outflow has essentially one 
degree of freedom besides the multiplicity M.. One can chose, e.g., p + as an independent 
variable, or any combination of p + and p_, ((p + ,p_), that is independent from Ai(p + ,p-). 
Below we will choose variable ( and formulate the dynamic equation for (; it will describe 
how the interaction with radiation governs the outflow dynamics along the magnetic field 
lines. 



5.2. Momentum equation 

Since the outflow is bound to move along the magnetic field lines, we will need the 
projection of momentum equation onto B. It is convenient to write this equation in a 
covariant form that is valid in any curved coordinate system x % (e.g. spherical coordinates). 
In a steady state, the momentum equation reads 

1^' = ^ < 2 °> 

where dP/dtdV is the rate of momentum exchange with radiation (per unit volume), is 
the covariant derivative (indices i, k run from 1 to 3), and T lk are the components of the 
plasma stress tensor, 

T k = m e c 2 I u*u k - = nm e c 2 ^ f V - f e (p) dp. (21) 
J 7 */ 7 



the slower negative charges still move freely. As p + deceases, the minimum-width model would require an 
increase in p_. Realistically, p_ would not first react to the reduced p + until particles with p = p_ also 
begin to experience the drag force (which is positive as p_ < p+). After this point, the outflow may not be 
far from the "minimum-width" state. 
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Fig. 8. — Relation between p + and p_ for the waterbag model shown in Figure 7. Each curve 
describes an outflow of a given multiplicity Ai, which is indicated next to the curve. 
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Here u l = pB l /B are the spatial components of the four- velocity vector of a particle with 
dimensionless momentum p, dn = nf e (p)dp is the number density of particles with momenta 
(p,p+dp), and dn/'y is the corresponding density in the frame moving with ^(p) = (1+p 2 ) 1 / 2 . 
Then the left-hand side of Equation (20) takes the form, 

§V fc r fc = ^V fc (nm e c 2 ^^ =m e c 2 B k V k [^), (22) 

where we used BiV \(B l B k / ' B) = (which follows from VkB h = 0). The directional deriva- 
tive B k Vk equals Bd/dl where / is length measured along the magnetic field line. Equa- 
tion (22) can be further simplified using n/3/BM = const, which follows from the relation 

n = M \ (23) 

ecp 

and j/B = const (cf. Equation 5). This gives 

Bi„ ^ _ , P d (M 



B v ' r -^ n jAji\T v n- (24) 

When Ai = const (i.e. no new pairs are created), Ai cancels out. Finally, the substitution 
of Equation (24) to Equation (20) gives the momentum equation in the following form, 

§ = J-„ C = f, (25) 
dl (3 m e c 2 (3 

where T = n~ x dPjdtdV is the average force exerted by radiation per particle. If T — 
(e.g. no interaction with radiation), ( remains constant along the field lines, which implies 
p± = const — the momentum distribution remains unchanged along the flow. 

For the waterbag model, (3 is given by Equation (17). The quantity p(3 can also be 
expressed in terms of p±, using the indefinite integral 



J p(3dp= j pd-y = y - i In (j~Z^J + const - 



This gives 



c = 



7+ -7- 



T 1 (l + /3 + )(l-/3_) 

o(P+7+ -P-7-) - T ln 



(26) 



(27) 



2^'^ 4 (1 -/5+)(l + 

In our numerical simulations, we use ( as the variable that describes the dynamical state of 
the outflow, and solve the differential Equation (25) for (. The values of p± that correspond 
to a given ( and M are found from Equations (19) and (27). 
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The force T experienced by an outflow with a given ( is determined by the local radiation 
field, which is described by the intensities in the two polarization modes I±(oj, n) and I\\(oj, n) 
(Appendix A). In general, the force acting on a plasma with a distribution function f e (p) is 
given by Equation (A20). Force T is easily calculated in the simple case of an optically thin 
magnetosphere exposed to the central blackbody radiation (Appendix B). Then the right- 
hand side of Equation (25) is a known function of ( and it is straightforward to numerically 
solve this equation; the results are similar to the two-fluid model described in Section 4. 
The optically thin approximation is, however, invalid for active magnetars. The radiation is 
not central; instead I± and I\\ must be calculated self-consistently, together with the outflow 
dynamics. This requires radiative transfer simulations. 

Note that we assume here that p± > 0, so that all particles move away from the star 
until they reach the equatorial plane and disappear (annihilate) there. This approximation 
is reasonable for the flow along the extended field lines with apex radii -R max ~ 10R. Near 
the apexes, i.e. near the equatorial plane 9 = n/2, the radiative drag enforces p± ~ p+ <C 1, 
creating a dense layer of slow particles where they annihilate (Fig. 2). A more general model 
of the flow could allow the particles to cross the equatorial plane and enter the opposite 
hemisphere. By symmetry, this would be equivalent to the reflection boundary condition, 
i.e. the mirror image of the outflow approaching the equatorial plane would emerge from 
the equatorial plane. Then the distribution function f e {p) must extend to negative p. This 
modification would be required for the plasma flow on field lines with small R max , where 
radiative drag is less efficient and the plasma can cross the equatorial plane with a large p. 
In this paper, we focus on the field lines with i? max ^ 10, where this does not happen. The 
field lines with small R max are assumed to form a cavity with j — and a negligible plasma 
density (Figure 1). 



6. Self-consistent radiative transfer 

The problem of radiative transfer in a relativistically moving e ± plasma whose velocity is 
controlled by the radiation field is not unique to magnetars. A similar situation may occur in 
accretion disk outflows (Beloborodov 1998) and gamma-ray bursts (Beloborodov 2011b). The 
strong radiative drag (measured by the coefficient T> 1, Equation [4]) was previously shown 
to simplify the problem, as it forces the plasma to keep the saturation momentum p* such 
that the net radiation flux vanishes in the plasma rest frame. This "equilibrium" transfer has 
one additional integral compared with the classical Chandrasekhar-Sobolev transfer problem 
for a medium at rest. The transfer in magnetar magnetospheres has, however, two special 
features that complicate the problem. First, the electric current and plasma instabilities 
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imply additional (electric) forces that broaden the momentum distribution around p*, as 
discussed in Sections 4 and 5. Therefore, the equilibrium condition p = p* is not satisfied 
even where Z> 3> 1. Second, opacity is dominated by resonant scattering, whose rate is 
sensitive to the particle momentum. Below we develop a method to solve the self-consistent 
transfer for magnetars. 

6.1. Monte-Carlo technique and the "virtual beam" method 

Radiative transfer in a magnetosphere filled with plasma with given parameters can be 
calculated using the standard Monte-Carlo technique (e.g. Fernandez & Thompson 2007; 
Nobili, Turolla & Zane 2008). Blackbody photons are injected into the magnetosphere at the 
neutron-star surface and their trajectories are followed until they escape the magnetosphere. 
We implement this method using the scattering opacity given in Appendix A and keeping 
track of the photon polarization, which can switch in the scattering events. 

Calculation of transfer in an outflow with self-consistent dynamics is a more ambitious 
goal, as the plasma parameters are not known in advance. A natural approach is iterative. 
One can start with a trial outflow, calculate radiative transfer to find the radiation intensity 
that would correspond to this outflow, and then re-calculate the outflow dynamics in the 
obtained radiation field. These iterations can be repeated until they converge. The problem 
is simplified for the waterbag plasma model (Section 5) as the outflow is described by one 
dynamic variable (. As the first trial one can take the outflow solution ((r, 6) for the optically 
thin magnetosphere exposed to the central blackbody radiation. 

One then encounters the following difficulty For the calculation of next iteration of 
outflow dynamics, one needs to know ij_(r, u>,n) and I\\(r,u,n) everywhere in the mag- 
netosphere. In axisymmetric magnetospheres, the radiation intensity is a function of five 
variables: two for location (radius r and polar angle 6), one for spectrum (u), and two for 
angular distribution (unit vector n is described by two angles). To determine intensities I± 
and J|| with sufficient accuracy, one has to introduce a five- dimensional grid and accumu- 
late large photon statistics for each grid cell during the Monte-Carlo simulation of radiative 
transfer. A grid of size N for each of the five variables has N 5 cells. Accumulation of photon 
statistics in each cell requires the calculation of a huge number of Monte-Carlo realizations 
of the photon trajectory in the magnetosphere. This is expensive and leads to poor accuracy. 

There is, however, a more efficient method that does not require the knowledge of in- 
tensities I± t \\(r,uj,n). They are not, in fact, needed for the calculation of outflow dynamics. 
What enters the momentum Equation (25) is the average force F(r, 9), which may be tab- 
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ulated in the two-dimensional space of r, 9. This force can be calculated directly during the 
Monte-Carlo simulation of radiative transfer if we find a way to evaluate the contribution of 
each simulated photon to F(r, 9) as it propagates through the magnetosphere. 

This can be achieved if we imagine that the photon is replaced by a beam of radiation. 
As we follow the photon along its trajectory, we can calculate the force that would be applied 
by the imaginary beam to any given outflow that can be imagined in the magnetosphere. 
In the waterbag model, the outflow is described by one dynamical parameter (, and so the 
force applied by the beam can be tabulated on a grid of (. 

Once we know how calculate the force created by each photon trajectory in our Monte- 
Carlo simulation, we should be able to average it over all simulated photons and thus ac- 
curately evaluate F(r, 9,C)- This gives the force that would be applied by our radiation 
field to an outflow with any given ( at any point r, 9. At the next iteration step T{r, 6, () 
is used to obtain the new outflow solution ((r,9) by integrating the momentum equation 
dC/dl = W(l,(), where W = F/fim e c 2 (Equation 25). Then the Monte-Carlo simulation 
can be repeated to calculate radiative transfer in the new outflow and find f"(r, 9, Q for the 
new radiation field. These steps can be repeated until the outflow solution £(r, 9) converges, 
i.e. remains practically unchanged by new iterations. 

The concrete implementation of this strategy is as follows. Let L th be the thermal lu- 
minosity of the star and N = L th /2.7kT be the number of photons emitted by the star per 
unit time. In our Monte-Carlo simulation we follow K MC ~ 10 7 random photon trajectories. 
This can be thought of as dividing N into Kmc random monochromatic beams. Each beam 
has a random start at the star surface (the photon energy is drawn from the Planck distribu- 
tion) and follows one random realization of the photon trajectory, which can involve multiple 
scattering events in the magnetosphere. The photon number flux in each monochromatic 
beam is N b = N/K MC , and the energy flux in the beam is 

E b = N b hw. (28) 

Note that N b = const along the beam (i.e. along the photon trajectory in the Monte-Carlo 
simulation) while the photon energy Huu changes after each scattering in the magnetosphere. 
The collection of K M c beams represent the state of the radiation field around the star. 

As we follow each realization of the photon trajectory, in parallel we calculate the force 
applied by the virtual beam to the plasma. We can imagine that the beam has a small cross 
section A (it will cancel in the final result) and flux density F = E b /A. Equation (A19) gives 
a general expression for the force exerted by radiation on a plasma with a given distribution 
function f e {p). For a monochromatic beam of frequency u propagating at angle $ with 
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respect to the magnetic field, this expression gives the following momentum deposition rate, 

2n 2 r e nF^[ ll f e (p 1 )- l2 f e (p 2 )], u S m$<u B , (29) 



dtdV e * u 2 

(dP/dtdV = if cusini? > w B ), where 




Pl , 2 (u>,#) = I cos^ T Jl - ^sin 2 ^ I , (30) 



uj sin 2 

and the factor £ depends on the beam polarization, 

1, JL 

1 ^ • 2 q (31) 
1 ^sm 2 $, h 

The rate of momentum deposition by the beam (i.e. the exerted force) per unit length along 
its trajectory is given by 

dP dP A 

■A. (32) 



dt ds dt dV 

We need to tabulate the force on a spatial grid, which is used to calculate the outflow 
dynamics at the next iteration. Therefore, we need to evaluate the net force applied to a 
given spatial cell. The number of particles in the cell is nV c where V c is the cell volume, and 
the force exerted by the beam per particle in the cell is given by 

f = ^^^«W.W-V.W1. (33) 

In the transfer calculation, we track the photon trajectory using small steps 5s, much smaller 
than the cells of the spatial grid. To obtain the force .F&(r, 9, () applied by the beam in a 
given cell we integrate Equation (33) along the photon path crossing the cell. Note that 
a finer spatial grid implies a larger dThjds oc V^T 1 ; however, it also implies that the cell is 
less frequently visited by photons in our Monte-Carlo simulation, and the photons spend a 
shorter time in the cell, therefore the final result does not depend on the grid. 

In an axisymmetric magnetosphere, the spatial cells are tori of volume Vij = 

2nrf sin 9j Ar A9. The net force applied per particle in a given cell is obtained by 

summing up the contributions Tb(ri,9jX) from all simulated beams, 

Hn, ^0 = -^-E^ / — * MM - Me( P 2)] ds. (34) 

K MC ^ V itj Jc C \\(i,j) U 

Note that the beam may cross a given cell multiple times in an opaque magnetosphere, and 
all crossings contribute to the path integral in Equation (34). The e ± distribution function 
feip) is determined by the parameter ( (assuming a given pair multiplicity Ai, see Section 5). 
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As a test, one can apply Equation (34) to the simplest case of a transparent outflow 
exposed to the central thermal radiation. Then the expected T can be directly calculated 
using Equation (B17) in Appendix B. We analytically verified that in this case Equation (34) 
is reduced to Equation (B17). We also tested our numerical code; it reproduced the analytical 
result. 



6.2. Results 

The obtained self-consistent solution for the outflow dynamics is shown in Figure 9. 
We show p + rather than our dynamical parameter (, because p + is closely related to the 
radiation emitted by the outflow. We find that only particles with the highest momenta 
p ~ p + resonantly scatter thermal photons in the relativistic zone p + > 1; the remaining, 
dominant part of the momentum distribution does not participate in scattering. The Lorentz 
factor of the scattering particles is 

7sc~(l+P+) 1/2 - (35) 

The solution p+(r, 8) shown in Figure 9 was calculated assuming that the star has radius 
R — 10 km, a uniform surface temperature kT = 0.3 keV, and a moderately twisted dipole 
magnetic field with -B po i e = 10 15 G and ip = 0.3; the multiplicity of the e ± flow was fixed at 
Ai = 200. The flow is injected at r = 2R with p + (2R) = 100. The choice of the boundary 
condition is not important as we are interested in the flow behavior outside ~ 5R, where the 
scattered photons avoid conversion to e ± pairs and can escape, i.e. where the observed hard 
X-rays are produced (see the accompanying paper B12). For any reasonable boundary value 
p + (2R), the flow relaxes to the same solution p + (r, 6) outside a few stellar radii. Remarkably, 
our simulations gave practically the same solution for a broad range of parameters tp, A4, T 
that is relevant for magnetars. This demonstrates a robust self-regulation mechanism; this 
important feature is explained below (see also the accompanying paper). 

There are two distinct zones in the outflow: 

I. Non-relativistic zone p + <C 1, which is near the equatorial plane (red zone in Figure 9). 
This zone has a huge optical depth and scatters essentially all thermal photons that impinge 
on it from the star (~ 10% of the thermal luminosity L th ). We will call this zone the 
"equatorial reflector." 

II. Relativistic zone p + > 1 (blue to green in Figure 9). This zone is transparent for 
essentially all thermal photons flowing from the star, except for rare photons in the far Wien 
tail of the Planck spectrum, which may be neglected. Basically, the outflow in this zone 
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Fig. 9. — Self-consistent outflow solution. Color shows the parameter p + of the waterbag 
distribution function. Hard X-ray emission discussed in the accompanying paper (B12) is 
produced by particles with p ps p + . The simulation assumed that the active j-bundle occupies 
the field lines with apex radii -R max > 10R (Section 1). 
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does not "see" the thermal radiation flowing from the star. It mainly interacts with (and is 
decelerated by) the quasi-thermal photons that flow from the equatorial equator. 

It is easy to see why the outflow in the relativistic zone 7+ 3> 1 is regulated so that 
it interacts with a tiny fraction of the thermal photons around the magnetar. Scattering 
on average boosts the energy of a thermal photon huu by a factor comparable to 7 2 , and 
hence the energy lost by the outflow per scattering is ~ rftvjj. If we imagine that each 
thermal photon is scattered once with an average blueshift of 7 2 , the generated hard X-ray 
luminosity 7 2 Lth would exceed the kinetic power of the outflow L, which is impossible. The 
scattering rate must be kept low, just sufficient for the self-consistent gradual deceleration 
of the outflow. The outflow "vision" of targets for scattering is controlled by the resonance 
condition 7(1 — f3 cos $)u = u B , where 1? is the angle between the photon direction and the 
particle velocity f3. The scattering rate is greatly reduced if 7 is so low that the resonance 
condition is satisfied only for rare photons in the Wien tail of the thermal spectrum. The 
situation may be compared with the regulation of the nuclear burning rate in the sun. The 
self-consistent temperature is sufficiently low so that fusion reactions occur only in the far tail 
of the Maxwell distribution; as a result only rare particles participate in burning, and these 
particles are concentrated in a narrow interval of (large) momenta, a phenomenon known 
as the Gamow peak. Similarly, our self-consistent outflow moves sufficiently slow so that 
resonant scattering is only enabled between rare particles with the highest p ~ p + and rare 
thermal photons with the highest energies in the particle rest frame. These lucky photons 
have the largest $ (gained after reflection from the equatorial reflector) and large fku 3> 
2.7kT. The inward direction (cos$ < 0) of the reflected photons gives them particularly 
large blueshift in the outflow rest frame and hence reduces the energy requirement in the lab 
frame. This makes the reflected photons the dominant targets for scattering even though 
their density is much smaller than the density of photons flowing directly from the star. The 
reflected photons have a diluted quasi-thermal spectrum of temperature T. 

In essence, we observe in our simulations that the outflow moves fast enough to reso- 
nantly interact with reflected photons of energy Tiw ~ (7— 10)kT (the low-density exponential 
tail of the thermal spectrum), and slow enough to not interact with the main peak of the 
reflected thermal spectrum ftio ~ 3kT. This condition, together with the resonance condition 
7(1 — f3cos , d)uj = ub and cos-# ~ —0.5, determines that the scattering plasma moves with 
Lorentz factor 

m e c 2 B 

7sc "TblTV (36) 

as long as 7 SC ^> 1. In this regime, the number of target photons visible to an outflowing 
particle has a strong exponential sensitivity to the particle momentum p. Essentially all 
scattering must be done by particles with the highest momenta in the distribution function, 
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i.e. with p ~ p+ for the waterbag distribution, as indeed observed in our numerical simula- 
tion. Therefore, 7 SC is associated with p + . Equation (36) serves as a simple and reasonably 
accurate approximation to the exact numerical results shown in Figure 9. Its applicability 
is not limited to the specific simulation with its -B po ie, T, vp, Ai — the approximation works 
well for other magnetar parameters, because of the robust self-regulation effect described 
above. This fact is further discussed and illustrated in Figure 2 in B12. 

Besides Equation (36), the transfer problem is characterized by the position of the 
equatorial reflector. It is described by a simple formula, which can be used to scale the 
results shown in Figure 9 to models with other parameters. The formula is based on the 
following fact (see Appendix B): if a given active magnetic loop extends to the region where 
fojjB < 20/cT, the central thermal radiation exerts a sufficiently strong drag on the outflow 
to bring it to rest at the top of the loop. The region ftw B < 20kT corresponds to r > Ri 
where 



Here ft = R 3 B po \ c /2 is its magnetic dipole moment of the star. Equation (37) describes the 
position of the inner edge of the non-relativistic (red) zone; e.g. the model in Figure 9 has 



It is instructive to compare the result of the full transfer calculation in Figure 9 with 
the simplest, optically thin two- fluid model in Figure 3. The relativistic zone p + > 1 remains 
practically transparent to thermal radiation. The key difference is the presence of the opaque 
equatorial reflector. The reflector weakly affects the spectrum of thermal photons supplied 
by the star, however it significantly changes their angular distribution in the magnetosphere. 
As a result, the radiation exerts a stronger drag on the outflow and p + decreases faster 
along the magnetic field lines. In Figure 3, 7 + ~ p + ^> 1 remains huge near the magnetic 
axis — the central radiation is unable to decelerate the plasma because the photons flow 
from behind and have small angles d with respect to the plasma velocity. In Figure 9, the 
equatorial reflector supplies photons with large i?, which efficiently decelerate the outflow, 
according to Equation (36). 

The upscattered photons of energy E ~ 7 2 £ t are beamed along the relativistic outflow. 
Therefore, they become unable to decelerate the plasma, even though they can scatter mul- 
tiple times before escaping. The outflow significantly loses energy when it scatters a photon 
propagating at a large angle d with respect to the outflow velocity; only in this case the 
scattering boosts the photon energy by the factor of ~ 7g C 1. After the scattering, the 
photon angle is reduced to •& ~ 7" 1 , and its subsequent scatterings have a small effect on 
the outflow dynamics. The beamed radiation initially moves together with the plasma and 




(37) 



Ri 



WR. 
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then escapes. Our transfer simulations include all scattering events, however practically the 
same p + (r, 6) would be obtained if only single scattering were allowed in the relativistic zone 
P+ > 1- 

7. Conclusions 

This paper examined the behavior of the relativistic plasma created by e ± discharge 
around magnetars. Plasma circulation in the magnetosphere is schematically shown in Fig- 
ure 2. We focused on large magnetic loops, which must be heavily loaded with e ± pairs. The 
plasma momentum is controlled by the radiation field around the star, which interacts with 
e + and via resonant scattering. We developed a method to calculate radiative transfer 
in the self-consistently moving plasma and obtained the solution for the e ± flow (Figure 9). 
The solution is a strong attractor — the behavior of the plasma outside a few stellar radii 
does not depend on how the plasma is injected near the star and its initial Lorentz factor 
7o, as long as 7 3> 30 (70 ^ 10 3 is expected, which corresponds to the discharge voltage 
$ e ~ 10 9 V). The e ± flow shows the following features: 

(1) The relativistic flow scatters radiation with a well defined Lorentz factor 7 SC given 
in Equation (36); 7 SC decreases proportionally to B along the magnetic field lines. 

(2) The relativistic flow remains transparent to thermal photons emitted by the star 
until it decelerates to non- relativistic momenta p < 1. The non-relativistic zone is opaque 
to the star radiation and forms the "equatorial reflector" (red zone in Figure 9). 

(3) The energy lost by the decelerating flow is converted to hard X-rays. The resulting 
emission is calculated and compared with observations in the accompanying paper (B12). 

(4) The plasma is nearly neutral, n + ps n_, and carries the electric current j = 
(c/47r)V x B by adjusting the particle velocities. In large magnetic loops (extending to 
the region of B ^ 10 13 G) the plasma has a high multiplicity Ai ~ 10 2 , and both elec- 
trons and positrons outflow from the neutron star, with a small separation in the velocity 
space. This separation is sustained by a modest electric field induced along the magnetic 
field lines. 

(5) The enforced electric current and radiative drag together create a configuration that 
is prone to two-stream instability, which is expected to generate low-frequency radiation. 
The mechanism of initiating the two-stream instability is unique to magnetars, explaining 
their special radio emission and possibly optical/UV emission. 
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A. Resonant scattering 

Resonant scattering plays a significant role in ordinary pulsars (e.g. Kardashev, Mitro- 
fanov, & Novikov 1984; Daugherty & Harding 1989; Sturner 1995; Lyubarsky & Petrova 
2000). It is also the dominant radiative process in magnetar magnetospheres, which governs 
the radiative transfer calculated in this paper. Below we summarize basics of resonant scat- 
tering, write down the cross section, the optical depth of the e ± flow, and the radiative drag 
force that are used in our numerical simulations. 

Photons scattered in the region B > 10 13 G are immediately absorbed (B12). Interesting 
radiative transfer occurs outside this region, where B <C Bq. In particular, the observed 
hard X-rays of energy up to several MeV are radiated where B <C Bq (B12). Electron recoil 
is small for resonant scattering in such relatively weak fields, and the scattering cross section 
is particularly simple. 



A.l. Scattering cross section 

In classical language, the electro- magnetic wave (photon) with frequency ujb = eB/m e c 
resonates with the Larmor rotation of electron. Then the wave strongly accelerates the 
charged particle and generates scattered radiation. The corresponding cross section is largest 
for waves with the right-hand circular polarization e = e_ that matches the electron Larmor 
rotation. Here e_ = 2 _1//2 (e x — ie y ) and e x , e y , e z is a Cartesian basis with the z-axis 
anti-parallel to B. For a wave with an arbitrary polarization vector e, only the projection 
of e on e_ is responsible for the resonance, and the cross section is reduced by the factor 
|e* -e„| 2 , where e* is the complex conjugate of e. In quantum language, the resonance occurs 
because the photon energy matches the energy Hojb needed for the electron transition from 
the ground Landau state to the first excited state. The resonance has a finite width. It 
equals the natural width of the cylcotron line T, which corresponds to the lifetime of the 
excited electron to spontaneous transition back to the ground Landau state. 

For an electron at rest, the differential cross section for photon scattering into solid 
angle dVt' is given by (Canuto et al. 1971; Ventura 1979), 

2 ^ I* |2 I /* |2 / A 1 \ 

m> =r * {uj-u B y + {T/2y |e ' e -' |e ' e -' ' (A1) 

where r e = e 2 /m e c 2 , uj is the photon frequency, e and e' are the polarization vectors of the 
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photon before and after scattering. Equation (Al) retains only the resonance peak of the 
cross section and neglects the non- resonant part. Positron cross section is described by the 
same equation except that e_ is replaced by e + = 2~ 1 / 2 (e, E + ie y ) (positrons gyrate in the 
opposite sense). The cross section can also be derived in the framework of quantum electro- 
dynamics (Herold 1979; Daugherty & Harding 1986). For Hlob *C m e c 2 (which corresponds 
to B <C Bq) the result is reduced to Equation (Al). 

The resonance line is very narrow, Y/lob = {A/Z)ol{B / Bq) <C 1, where a = e 2 /hc = 
1/137 (Daugherty & Ventura 1978; Herold et al. 1982), and the resonance factor \{uj — ujb) 2 + 
(r/2) 2 ] -1 is well approximated by the delta- function 2tiY~ 1 5{uj — ujb)- Then the cross section 
may be written as 

% = 2tt ^ = 3n 2 r e c5(uj - u B ) |e* • e ± | 2 |e" • e ± | 2 , (A2) 
a/r ail' 

where +/— correspond to scattering by positron/electron, fi' = cosi?', and i?' is the angle 
of the scattered photon with respect to the magnetic field B. The distribution of scattered 
photons is axially symmetric about B; this fact has been used in Equation (A2). 

Two polarization states (eigen modes) exist for the photon. They are controlled by the 
dielectric tensor of the magnetosphere. In the considered region r < 100-R, the dielectric 
tensor for photons of interest (X-rays) is dominated by the magnetic vacuum polarization 
effect (e.g. Beresteskii et al. 1982); the plasma contribution to the dielectric tensor is much 
smaller and may be neglected. Magnetic vacuum defines two linearly polarized eigen modes 
for electromagnetic waves: which is perpendicular to the (k, B) plane, and ey which is 
parallel to the (k, B) plane (here k is the photon wave vector and e shows the direction 
of the electric field in the wave). The _L and || modes are also called E-mode and O-mode, 
respectively. Their refraction indices are 

^= i+ ^M) 2 ^' n ^ i+ ^{^^ <A3) 

where d is the photon angle with respect to B. The two modes have slightly different 
propagation speeds c/N and therefore they adiabatically track, i.e. the photon propagating 
through the curved magnetic field preserves its polarization state. The adiabaticity condition 
reads kl B {N\\ — N±) 3> 1 where l B ~ r is the characteristic scale of the spatial variation of 
B (see e.g. Fernandez & Davis 2011 for a detailed discussion). This condition is satisfied 
for X-rays in the considered region of the magnetosphere where scattering occurs. Thus, in 
our transfer problem, the photon can switch its polarization state only in a scattering event. 

As the photon can be in either polarization state, calculation of radiative transfer in- 
volves four scattering processes _L— >■_!_, _L— ||— >■_!_, and ||— The corresponding cross 
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sections are given by Equation (A2) with |e^_ ■ e ± | 2 = |e'j_* • e ± | 2 = 1/2, |ejj • e ± | 2 = yU 2 /2, 
and |ej|* • e ± | 2 = /i' 2 /2. Note that the cross sections of electron and positron are equal, as 
|e* • e_| 2 = |e* ■ e + | 2 for any linear polarization e. 

Equation (A2) describes the cross section of electron (or positron) at rest. In our transfer 
problem, the particles are moving along B and the above equations should be used in the 
rest frame of the particle. Note that the polarization states and ojb are invariant under 
Lorentz boosts along B. Consider a photon with energy hu and propagation angle with 
respect to B, in the _L or || polarization state. We are interested in its scattering by an 
electron (or positron) that moves with velocity (5c along the magnetic field line. The total 
cross section in the lab frame may be obtained by integrating the differential cross section 
in the electron rest frame and then multiplying the result by 1 — /3cos$, 

a tot = 27T 2 r e c^5{u -u B )(l- Pfi), (A4) 

where 

u> = 7 (1 (A5) 

is the photon frequency in the electron rest frame and fi = cos$. The factor £ depends on 
the photon polarization and is given by 

where \x = cos'd and d is the photon angle with respect to B in the electron rest frame. The 
cross section <7 to t is summed over the final polarization states. 

The outcome of the scattering may be either _L or || photon propagating at angle d' in 
the electron frame. The probability distribution for ji' = cos$' is given by 

W = ££ = !r. (AT) 

where £' = 1 or jl' 2 , depending on the final polarization state. The integral / P(fi') dp,' 
equals 3/4 for the _L state and 1/4 for the || state. Thus, 3/4 of scattering events produce _L 
photons with uniform distribution P{fa') = const and the remaining 1/4 of scattering events 
produce || photons with P(jl') oc jl' 2 . The distribution of the final photons over angle and 
polarization does not depend on the initial state of the photon before scattering. 

The standard description of resonant scattering summarized above assumes the transi- 
tion between the ground and first excited Landau levels, which has the largest cross section. 
Transitions to higher levels are neglected in our transfer calculations. 
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A. 2. Opacity 

Optical depth of a relativistic plasma to resonant scattering was discussed previously in 
detail (e.g. Fernandez & Thompson 2007 and refs. therein). Here we write down relevant 
equations and introduce notation that is used below in the discussion of radiative drag. 
Consider a photon of energy froj that propagates through e ± plasma with density n = n + +n_. 
The plasma particles move along B with momentum distribution f e (p) that is normalized 
to unity 

J f e ( P ) d P = i. (A8) 

The optical depth dr seen by the photon as it propagates distance ds is given by 

dr f 

fa= n °"t°t /e(p) d P- ( A9 ) 

Assuming that the magnetosphere has a small optical depth to non-resonant scattering 
(Appendix C), we include only resonant scattering and use Equation (A4) for er tot . The 
integration over p may be carried out using the identity, 

where 

^ = 5(7-w)« = oj-mK (AH) 

and pk are all possible solutions of equation Cj{p) = ujb- The delta-functions S(p—pk) express 
the fact that the photon is scattered by electrons or positrons with momenta pt for which 
the resonance condition is met. The relation a)sini? = wsin-^ together with the resonance 
condition u = ub determines sm$ = (uj/ujb) sin?? and leaves two possibilities for ft, = cos$, 



J 2 \ 1/2 

fti = ± ( 1- — sin 2 ^ . (A12) 



Angle d exists (i.e. the resonance is in principle possible) for photons that satisfy the 
condition ousmd < ug. Then Equation (A12) defines two electron velocities flip, which may 
be found from the Doppler transformation of the photon angle, p, — (fx — j3) / (1 — (3fi). It 
yields, 

P = f^4, A,2 = f^. (A13) 
1 — /jifj, 1 =F WW 

The corresponding Lorentz factors 7 = (1 — /3 2 )~ 1/2 and dimensionless momenta p = 7/? are 

71,2 = . \*r~ , Pi,2 = f^ . (A14) 
sin v sin v sin v sin v 
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Substitution of Equations (A4) and (A10) to Equation (A9) and integration over p gives 

^ = 2tt 2 r e - -1 n if e ( Pl ) + fM] • (A15) 
as ui 

Here £ = 1 for _L photons and C, = fl 2 for || photons (Equation A6). 



A. 3. Radiative drag force 

Consider now a radiation field with intensity I(u, k) in a given polarization state, _L or 
|| . As a result of scattering, radiation exerts a force on the e ± plasma. We are interested in 
the component of this force along the magnetic field. The force applied to unit volume of 
the plasma is given by 



dtdV 



Jdnjdujdp n f e (p) a tot AP. (A16) 



Here J dVt is the solid-angle integration over photon directions n = k/fc, and AP is the 
average momentum (per scattering) passed to an electron or positron with Lorentz factor 7 
by a photon (cu,k). In the electron rest frame, AP equals the photon momentum along B, 
as its average momentum after scattering vanishes (resonant scattering is symmetric in the 
electron frame when huj B <C m e c 2 ). Thus, AP = jifuJj/c and the Lorentz transformation of 
the four-momentum vector to the lab frame gives 

AP = 7 ^, (A17) 

c 

where we used the condition oj = ub since we consider only resonant scattering. 

Radiation is described by two intensities I±(oj,\s) and I\\(u,k) in the two polarization 
states. They scatter with cross sections <7 to t that differ by the factor jl 2 (see Equations. A4 
and A6). Substituting Equation (A4) to Equation (A16) and taking the sum over polariza- 
tions, one finds the net force exerted by I± and Jy, 

dP f Jr ^ f J f J (I ± + fi 2 I\ 



J dQ J du J dp ^ J " - nf e (p) 2n 2 r e c5{u - u B ){l - pp) AP. (A18) 



dtdV 

Integration over dp similar to that in Section A. 2 gives 

dP r ruiBsintf 



dtdV 



/ i>u>b smtt 
dn dco 27t 2 r e -| (7 ± + /i 2 /||) n [ 7 i/ e (pi) - 7 2 /e(p 2 )] , (A19) 



where pi )2 (c<;,^) and 7li2 are given by Equations (A14) and (A12). The upper limit in the 
integral over u takes into account that only photons with u < ujb^^ m &y be resonantly 
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scattered (Section A. 2). Equation (A19) is useful because it shows the contribution of each 
photon (u, k) (in the _L or || polarization state) to the drag force. It can be used even if the 
radiation intensity is not known in advance and needs to be found self-consistently with the 
plasma dynamics (Section 6). 

An alternative way of simplifying Equation (A18) is to first integrate over uj, which gives 

^ = JdnJ dp27r 2 r e (I ± + jl 2 I ll )nf e (p) 7(^-/3), (A20) 

where I± and I\\ are evaluated at uj — 7 _1 (1 — f3fi)~ 1 u B . Equation (A20) is convenient to 
use when the radiation intensity is known. 



B. Optically thin outflow 

Consider a magnetosphere that is optically thin to resonant scattering, so that intensity / 
is dominated by the unscattered radiation from the star. We will assume that the neutron star 
emits approximately blackbody radiation with the _L polarization (the _L photons dominate 
the surface radiation because they have a larger free path below the surface, e.g. Silant'ev 
& Iakovlev 1980). Then 

^ h = 0, (Bl) 



87i 3 c 2 [exp(fkj/kT) - 1]' 



and we deal with the outflow dynamics in the known radiation field. This case was studied 
in detail in previous work (e.g. Sturner 1995). 



B.l. Scattering rate for one particle 

Consider an electron (or positron) located at r, 9 and moving outward with Lorentz 
factor 7 = (1 — /3 2 ) -1 / 2 along a magnetic field line. The number of photons scattered by the 
electron per unit time is 

Nx = [ dtt f du'-^la^ /( "" [nl ' n) dti, (B2) 

J J nw 7 J riw res 

where we substituted cr to t (eq. A4). The resonant frequency depends on the photon direction, 
u TCS = 7~ 1 (1 — (3 •■n)~ l uj B where n = k/fc is the unit vector corresponding to dfl. 

At large radii r ^> R all photons at a given location r have approximately the same 
direction n || r. The angle between the stellar photons and the particle velocity, is given 
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by 11 — cos-# ps S r /B (assuming i? > -R/r). Then and integration over c?0 in Equation (B2) 
is reduced to multiplication by AQ ps n(R/r) 2 , the solid angle subtended by the star when 
viewed from radius r. This gives, 

Nsc= 2 ^[pA Urcs = (B3) 
x 2 7 ftw rcs 7(1 - p» 

where x = r/R. It is instructive to write iV sc in the following form, 

sc 4x 2 7 A y ' 1 j 

where 6 = kT/m e c 2 ps 10~ 3 and g(y) is the dimensionless Planck function evaluated at the 
frequency w rcs = 7 ~ 1 (1 - 

= * = It = 7 (i-/W (B5) 

where b = B/Bq. 



B.2. Drag force exerted on one particle 

The drag force applied by the central blackbody radiation to the electron is J 7 = N SC AP 
where AP is given by Equation (A17). It yields, 

^(/3) = ^^e 3 7 ^)(A-/3), (B6) 

where /3* = //. Force J 7 vanishes if /3 — /3+; in this case the radiation flux measured in the 
rest frame of the particle is perpendicular to B and cannot accelerate or decelerate it. In 
a weakly twisted magnetosphere, the magnetic field in the outer corona is approximately 
dipole. Then the saturation velocity /3* at a point r, 9 (spherical coordinates) depends only 
on 6 and is given by 

n B r 2 cos 6* 2 cos 6* . . 

B (1 + 3 cos z 6) L/ ' L suit/ 

The radiative force always pushes the particle toward p = p+. This effect may be measured 
by the "drag coefficient," 

Z> = f. B8 

c p dt v 7 
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Consider an electron (or positron) with momentum p ps p± = 7*/3*. A small deviation p — p± 
causes drag V oc p — p±, which may be written as 

V=vJl-2-\ (B9) 



where 

4 r e x 7^ \0.5 keV / 

Here = 67^/6 corresponds to photons that are resonantly scattered by the electron with 
p ps p*. The momentum is a strong attractor if X>* ^> 1. The value of X>* is sensitive to 
y+. In particular, in the equatorial plane, we have 7* = 1 and 

-,4 ( -Bpoic "N / kT \ ( r 



fcMl . 6xI0 .^^__J , (, = lr/2) , (Bll) 

where -B po i e is the dipole field at the magnetic pole. The condition V* > 1 corresponds 
to y± ^ 20. This implies that the e ± flow on magnetic field lines extending sufficiently far 
from the star is stopped by the radiative drag in the equatorial plane. For typical magnetar 
parameters, the flow stops on field lines with -R max ~ 10R. 



B.3. Optical depth in the single-fluid approximation 

The single-fluid flow has a distribution function f e {p') = 5(p' — p) where p(r,9) is the 
flow momentum. Then any photon of energy fvjj may only be scattered on the infinitesimally 
thin resonance surface defined by 7(1 — (5p)u = uub, where \x = cos^ describes the photon 
angle relative to the flow velocity. The optical depth of the resonant surface may be obtained 
from Equation (A9), which gives 

dr 

— = 2n 2 r e cnZ(l-f3fi)5(u-u B ), (B12) 
where u = 7(1 — f3fi)u. One can use the identity, 

*(*-**) = E 1 ff!~ 8fc \ r (B13) 

The location s& on the photon trajectory is where the photon crosses the resonant surface. 
Performing the integration over s along the photon trajectory, one finds the optical depth 
for one crossing of the resonant surface, 

27r 2 r e cn£(l - f3p) 



{u - u B )\ 



(B14) 
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If we specialize to the case of central photons emitted by the neutron star with the _L 
polarization, 6 then £ = 1 and fi = B r /B. For a moderately twisted dipole magnetosphere, 
the electric current density is given by j cip _B/47ri? max (Beloborodov 2009), and 

M ev 4vre/3 J R max ' ( ' B15 ' ) 

Note that n is small on field lines with a large -R max , which implies a low optical depth near 
the axis; this fact was first emphasized by Thompson et al. (2002) who used a self-similar 
twist model. 

The expression for the optical depth becomes particularly simple if the outflow has 
the equilibrium momentum p = p±. Then du/ds = 0, i.e. u remains constant along the 
radial ray through the outflow. This fact can be derived by noting that the Doppler factor 
7(1 — Pfj,) = 7" 1 is a function of 9 only — it does not depend on s = r for the approximately 
dipole magnetosphere. It is also easy to see that duiB/ds = — 3u>b/t, and we find 

7T lA , sin 4 0(1 + 3 cos 2 e) 1 ' 2 , 

The single-fluid model with p = p± may approximate the outflow only sufficiently close 
to the equatorial plane where 1 — /?* ^ M^ 1 . Nevertheless, the approximate Equation (B16) 
shows a general feature: the optical depth seen by the central photons is dramatically in- 
creased toward the equatorial plane (r oc sec#) and dramatically reduced toward the axis 
(r oc sin 9). Thus, a distant observer can see the unscattered radiation from the neutron- 
star surface when the line of sight is within a moderate angle 9 < tt/4 from the polar axis. 
This feature becomes even more pronounced in the full radiative transfer problem where the 
relativistic outflow is decelerated by the reflected radiation from the outer corona. Then 
scattering of the central radiation is negligible in the entire relativistic zone of the outflow. 



B.4. Drag force exerted on a plasma with a broad distribution function 

Equation (B6) describes the drag force exerted by the central thermal radiation on a 
particle with a given momentum p = 7/?. One can also consider a collection of particles with 



6 The neutron-star radiation is dominated by the _L polarization (Silant'ev & Iakovlev 1980). In addition, 
for the outflow with the equilibrium momentum considered below, the scattering of the || photons (even 
if they were included) would be suppressed. In the outflow rest frame, the photons move perpendicular to 
B, and the resonant cross section for the || polarization mode vanishes. 
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a momentum distribution f e (p) and derive the average force per particle T = n 1 (dP/dVdt). 
Equation (A20) gives 



where y is given in Equation (B5). The same result is obtained by averaging the force J 7 
given by Equation (B6), T — J J r (p)f e {p) dp. 

Equation (B17) simplifies when the plasma is described by the waterbag distribution 
function f e (Section 5.1); it leads to a straighforward calculation of the flow dynamics in 
the central radiation field. We use this simple outflow model as the first trial to initiate the 
iterations that converge to the solution shown in Figure 9. In the final solution, the drag 
exerted by the central radiation turns out negligible in the relativistic zone; instead, the 
outflow deceleration is controlled by the radiation streaming from the equatorial reflector, 
as discussed in Section 6. Then the force T derived in this section may be of interest only 
in the non-relativistic zone. 



Non-resonant scattering is not limited by the resonance condition, and hence many 
more photons can participate in scattering, although with a smaller cross section. Below we 
discuss the effect of non-resonant scattering on the dynamics of e ± flow around magnetars. 

Non-resonant scattering occurs mainly with photons in the Wien peak of the ther- 
mal radiation flowing directly from the neutron star, which dominates the photon density 
around the star. Relativistic particles see the thermal photons (of typical energy E ~ 3kT) 
blueshifted as E = 7(1 — f3p)E where p = cosfi describes the photon direction relative to the 
particle velocity in the lab frame. Sufficiently far from the star (where R 2 /r 2 < 1 — B r /B) 
the radiation can be approximated as a narrow radial beam; then p = B r /B. In general, 
p is a function of the particle position r, 9, <fi in the magnetosphere. For an approximately 
dipole field, p = 2 cos# (1 + 3 cos 2 9)~ 1 / 2 is a function of the polar angle 9 only. 

Magnetic field strongly affects the non-resonant scattering cross section if E < Hoob = 
bm e c 2 . If the electron is relativistic, the target photons are aberrated in the electron rest 
frame, cos-$ = p = (p — j3)/ (1 — flp). In the limit /3 — > 1, even photons with the || polarization 
have electric fields almost perpendicular to B, which makes their scattering inefficient. For 
photons with E Jtwb, the non-resonant scattering cross section is given by (e.g. Canuto 
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C. Non-resonant scattering 
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et al. 1971) 

a\\ I E \ sm z d a± IE 




E<^hu B . (CI) 



We assume <C m e c 2 and neglect Klein-Nishina corrections. Most of the radiation emitted 
by the neutron star has the _L polarization. 

The energy loss of the electron due to scattering is given by 

E e = - J dtt J dE(l-f3 f i)a(E) 1 -^^- (W-E), (C2) 

where n is the unit vector describing the photon direction in solid angle dft, and E' = jE is 
the mean expectation for the photon energy after scattering. This gives, 



E e = - J dQ J dE(l- /3/i) [ 7 2 (1 - /3/i) - 1] J(n, E) a(E) dE. 



(C3) 



In the simplest case of Thomson scattering of isotropic radiation, averaging over random \x 
gives the standard result E e = — (4/3)<tt c U 7 2 /3 2 , where U is the energy density of radiation. 
In our case, a < <7t, and the radiation field is not isotropic; far from the star it is better 
approximated central beam. 

Using Equation (C3) one can show that non-resonant scattering makes a small contribu- 
tion to the radiative drag compared with resonant scattering, and hence its inclusion in the 
calculation does not significantly change the outflow solution shown in Figure 9. Consider 
first the non-relativistic zone p + < 1. An upper bound on the non-resonant E e is obtained 
if we substitute in Equation (C3) cr(E) = <t t and /i — 0. This gives, 

E e = a T cp 2 U. (C4) 

The drag coefficient due to non-resonant scattering is defined similar to Equation (B8). Using 
U ~ L t h/47rr 2 c and dp/dt = E e /(3m e c 2 , one obtains 

V » / TLth7 „ » 0.2 L th , 35 r^ 1 7 < 1. (C5) 
47T r m e & 

In the relativistic zone p + > 1, the upper bound given by Equation (C5) increases propor- 
tionally to 7 and becomes useless, because it does not take into account the strong reduction 
of the scattering cross section below or- In this zone, the outflow adjusts so that it can 
resonantly scatter photons with E ~ 7kT and fi ~ —0.5 (photons flowing from the equa- 
torial reflector). This implies that the main targets for non-resonant scattering (photons 
flowing from the star with E ~ 3kT and /i > 0) have energies well below the resonance 
energy, E ~ (0.1 — 0.2)fruiB- Then the scattering cross section is strongly reduced below a T 
according to Equation (CI). When this reduction is taken into account, one obtains V < 1. 
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